Phase diagrams of a p-wave superconductor inside a mesoscopic disc-shaped sample 
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We study the finite-size and boundary eff'ects on a time-reversal-symmetry breaking p-wave su- 
perconducting state in a mesoscopic disc geometry using Ginzburg-Landau theory. We show that, 
for a large parameter range, the system exhibits multiple phase transitions. The superconducting 
transition from the normal state can also be reentrant as a function of external magnetic field. 
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Studies of Fermi superfluids and superconductors with 
multi-component order parameters have drawn much at- 
tention over the last few decades. A well-established ex- 
ample is spin-triplet p-wave superfluid '^He 1], the order 
parameter of which is a complex 3x3 matrix instead of a 
single component scalar as in the case of s-wave supercon- 
ductors. Many unusual behaviors of this superfluid are 
known, in particular "textures", where the order param- 
eter varies in space in a non-trivial way due to external 
fields, flows, or confining walls. There have also been 
intense studies of superconductors that are believed also 
to possess multi-component order parameters, for exam- 
ple, UPts and Sr2Ru04 [3]. While the precise order 
parameters in both cases are still controversial, many 
believe that both these two superconductor have order 
parameters which break time-reversal symmetry. Bro- 
ken time reversal symmetry necessarily requires a multi- 
component order parameter. A single component order 
parameter, though it can belong to a non-trivial one- 
dimensional representation, cannot break time-reversal 
symmetry in the sense that its complex conjugate differs 
from the original one only via a gauge transformation. A 
superconductor with broken time reversal symmetry can 
have exotic properties, such as circular dichroism and 
birefringence (Kerr rotation) [J| , internal magnetic fields 
[5| and surface currents 6] , to name a few. Experiments 
claimed to support broken time reversal symmetry have 
been reported both for UPta 0,3 and Sr2Ru04j9, 
though negative results are also in the literature [111 . 
Other aspects of intense recent interest are half-quantum 
vortices [IjI arising from the spin degrees of freedom, and 
Majorana vortex bound states [1J|, which would be pos- 
sible if the order parameter is a p-wave with the form 
Px + ipy, which is the case proposed for Sr2Ru04. 

In this paper, we study a two-component superconduc- 
tor in a confined geometry. Our motivations are several 
folded. First, as mentioned, in the context of '^He, a 
confining geometry can induce a non-trivial texture. A 
surface is necessarily a strong breaker of rotational invari- 
ance, and hence its effect on the order parameter depends 
on the relative orientation between the two. The ener- 
getically most favorable configuration therefore does not 
necessarily correspond to simply taking the uniform bulk 
order parameter and suppressing its magnitude near the 
surface. Second, for a multi-component order parame- 
ter, or more precisely an order parameter that belongs to 



a multi-dimensional representation, the different compo- 
nents possess the same transition temperature T^ in the 
bulk in the absence of external perturbations (by defini- 
tion). However, external perturbations can split the de- 
generacies, resulting in multiple phase transitions. This 
has been discussed in the context of thin-films of ■^He- 
B [15], as well as for UPta under the influence of the 
underlying anti- ferromagnetic order [2]. There have also 
been many recent experimental studies of mesoscopic su- 
perconductors, s [1^1, p [13], and d wave [l3|, includ- 
ing some interesting theoretical predications for the lat- 
ter, e.g. [13, but the physics associated with the multi- 
component nature of the order parameter is less explored 
in the literature. 

We shall consider a thin circular disc of radius R lying 
in the x-y plane. Variation of the order parameter along 
z, as well as the magnetic field generated by the super- 
current of the sample, will be ignored. We shall study a 
superconductor with a two-component order parameter. 
The two components r]x,y are supposed to transform as 
a vector under rotations within the x-y plane. We shall 
study how the order parameter varies over the disc. For 
definiteness, rj^.y are taken to be the two in-plane compo- 
nents of the orbital part of the p-wave order parameter, 
that is, the momentum dependence is rjxPx '^VyPy How- 
ever, we expect that many of our findings should be com- 
mon to other superconductors with multi-component or- 
der parameter. T his p oint will be discussed again below. 
Recently, a group '20] has studied theoretically this same 
system using Bogoliubov-deGennes equations. Their re- 
sults however differ significantly from ours. A comparison 
will be given later. 

We shall employ Ginzburg-Landau (GL) theory, but 
as we shall argue later, our conclusions are more general. 
The GL free energy density (per unit area) T consists of 
several contributions. The bulk contribution, J-'h, can be 
written as 



Jfc = a{T]* ■ rf) + I3i{r]* ■ rf) 
where a — a'{t — 1) with a' > 0, t 



Vp2\ri-vV, (1) 

T/T^ is the ratio of 
the temperature T relative to the bulk transition temper- 
ature T", 77 = rixX -\- rjyij, which we would often denoted 
as {rix,rjy). Stability requires pi > 0, /3i > —(i2- We shall 
take /32 > 0, so that for the bulk the equilibrium order 
parameter below T° has the form fj — ^(1, ±i), so that it 
has broken time-reversal symmetry, with \^\ = ji^-Y^"^- 



In the presence of gradients, there is an additional con- 
tribution to J- given by 

T, = K,{D,m){Djf^ir+K2{D,r^,){Dif^ir+K,{D,r^i){Dnjj 

(2) 
where Dj = dj + ^Aj is the gauge invariant deriva- 
tive. A is the vector potential, and the electron charge is 
— e. Repeated indices j, I are summed over x, y in eq.([2|). 
In writing down eq.([T]) and ^, we have ignored crystal 
anisotropics for simplicity, but we do not expect signif- 
icant qualitative change in our predictions below. The 
more general forms can be found in, e.g., [2l],[23. Stabil- 
ity requires Ki > 0, Ki23, = K1+K2+KS > 0. If we take 
rjx^y to represent the two in-plane orbital components of a 
pure p-wave order parameter, assume that the Fermi sur- 
face is isotropic in the plane, then, within weak-coupling 
theory, 132/ l3i = 1/2, and Ki = K2 = ^Ithe later holds 
up to particle- hole symmetric terms) [l|, Q , but we shall 
treat these coefficients as general parameters. 

We shall limit ourselves to solutions which are cylindri- 
cally symmetric, up to an overall gauge transformation. 
To this end, it is convenient to introduce the cylindrical 
coordinates (r, (/>) for space where (j) is the angle between 
r and x, and define r]± = rj^ i irjy. The bulk mini- 
mum energy solutions thus have ry_ = 0, r/+ ^ 0, or 
vice versa. r]± can be expanded as 77+ — J2n C+ {f')e™'^ , 

rj^ = X^„ C_ (r)e*'^"~^-''^, where n is integer and an ex- 
tra —2 is introduced in the rj- formula for convenience 

below. We have rf-f = i X;„(C+ ^ W + C'-"\^))e*'""^^'^ 

and7/-0 = -iiX;„(Ci"HO-Ci"H'^))e'^""^^'^- For solu- 
tions that are cylindrical symmetric up to a gauge trans- 
formation, only Cj_" for one particular n can be finite. 
Therefore different solutions are classified by n. In these 
cases, Cj_ can be chosen real without loss in generality. 
We shall first assume that the surface at r = i? is 
smooth. If the order parameter fj represents the mo- 
mentum part of the pairing wavefunction, the component 
perpendicular to the surface should vanish [23|, thus 



if ■ r — (r ^ R), 



that is. 



^(«) 



Ci"' - 



(r^R). 



(3) 



(4) 



The parallel component 7711 = rj ■ (f> should have vanishing 
gradient perpendicular to a planar surface [23|. At a 
surface with a finite curvature, it should satisfy [23| 



K-\-—m = — - 



dr 



K3 

r' 



R), 



(5) 



hence 
or 



Ki-M 



c 



in)-, 



K3 

R 



^(n) ^(n) 



(C|"^-Ci"0 (r^R)- (6) 



The differential equations satisfied by Cj_" (r) can be 
found by simple variation, noting that the total free en- 
ergy is F = 27rL drr{Tb+J^g)- We remark here that the 
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FIG. 1: Zero field Phase diagram: R is in unit of ^ = 
y^Ki23/a' ■ The parameters are /3 = /?2//?i — 0.25, ki = 
1/3 4- 5k, K2 = kz = 1/3 — 5k/2, with Ki = Ki/K\2s- Inset: 
R-t phase diagrams for different (5k's. 



boundary conditions eq.([3]|6|) guarantee that there are no 

net surface terms proportional to the variation 5Cj^ (R) 
of the order parameters at the surface, when we integrate 
by parts the gradient terms. They also guarantee that 
the normal component of the current vanishes for arbi- 
trary choice of order parameter profiles [M] , as should be 
the case of impenetrable walls at r = i?. 

First consider zero external magnetic field. It is worth 
having some analytic solutions before we show the nu- 
merical results. The transition temperature tj; from the 
normal to n = 1 state can be found analytically. The 
n = 1 solutions are time-reversal symmetric, and the 

free energy is invariant under C\_ < — > ±C_ . Since the 
boundary conditions are also symmetric under this trans- 



-^(^h^-\ = i^r'(i) 



.(1)n 



formation, the equations for Cs'' (r) = ^(C_j_ +C'_ ) and 
C], '(r) = ^{C-^ — Cl^) decouple. One possible solution 



is Ci^\r) = (hence ^- f = for all r), dj^\r) ^ 0, 

I.e. d+\r) = -Ci^V) = C^^Hr)- [25] After lin- 
earizing in the order parameter, it can be shown that 



the Bessel function of the 



C(i)(r) ex M^^illl^r), 

first kind. With eq.([6]), we find the relation 

Ki 2 



(l-t«) = 



(7) 



where a should satisfy aJ[{a) — ^Ji(a). The critical 
temperature is suppressed by a factor oc 1/R^. There- 
fore, we can find the critical radius via R^ — ay^Ki/a', 
which is the minimum radius of the system to maintain 
the superconductivity with n = 1 at zero temperature. 
Note that, due to eq.([7]), it is more convenient to set 
the vertical axes of R-t phase diagram to be l/R? , as 
shown in FiglTJ For this solution, we note that a — > 
as Kj, -^ Ki, which means that in the weak-coupling 
limit the transition temperature for our disc with smooth 
boundary is not at all affected by the finite radius and in 
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FIG. 2: Order parameters (a,c) of the ground states and the 
corresponding current distributions (b,d) for zero (black line) 
and finite field (red dash). 5k = 0, ^S = 0.25, i? = 6. The 
upper plots are for n = and i = 0; the lower plots are for 
n = \ and f = 0.7. The order parameters are normalized 
to ^Jol' / P\, the external field B is normalized as 6 = ^sS^^ 
(roughly to the bulk upper critical field), 5k and /3 as defined 
in the caption of Fig.l. 



fact is the same as that of the bulk, hence 1/i?^ ~^ o°- 
This can also be shown within the quasi-classical approx- 
imation, [26i| , and is therefore not an artifact of the GL 
approximation. The second case with analytic solution 
is K2 = K3 = 0. The order parameters near the critical 
— " -Jn{x) and C_ = Cn-Jn-2{x). Here 



point is C\ 



'-n+<- 



"'(!-*) 



^ . Using the boundary conditions ^ and 

(|n]), we can find the R-t relations for all n's. We find 
that, in this case, the state with the highest transition 
temperature is n = 0, which is degenerate with n = 2. 

Our obtained phase diagrams are summarized in FigU] 
Only n = 1 and n = (or equivalently, its time-reversed 

partner n = 2 with CS^ {r) -^ C'lp (^)) can be ground 
states. The left inset shows the kind of R-t phase dia- 
gram for < (5k < 0.44 27] . When the radius is smaller 

than a critical value. Re (y- intercept of red circle line) , 
the system does not have any superconducting phase. 
If the radius is between Re and R\^* (y-intercept of 
blue square line), which corresponds to the first order 
phase transition from n = 1 to n = with lowering 
temperature, it just has the superconducting phase with 
n = 1 after the second order phase transition (red cir- 
cle line). For the radius larger than i?^**, we have sec- 
ond order phase transition from the normal state and 
then a first order phase transition (blue square line) to a 
spontaneously time-reversal-symmetric-broken supercon- 
ducting state [n — Q). In the plot, we still show where 
the normal state would have become unstable toward the 
n = state by a dashed line, which does not correspond 
to a real phase transition for the system since it is below 
the n = 1 transition. On the other hand, if 5k is large 
enough (> 0.44) , the R-t phase diagram should be sim- 
ilar to the right inset of FigHJ The system just has the 
possibility to be in the superconducting state with n = 



at low temperature (the unphysical n = 1 instability line 
from the normal state is also shown as dashed). Transi- 
tion lines corresponding to second order phase transition 
from the normal state in R-t phase diagrams are linear 
within GL theory. For the first order transition lines, we 
found numerically that they are still practically linear. 
The main FigU] displays the critical radii at zero temper- 
ature. This figure can also be regarded as a plot of the 
critical radii at finite temperatures after rescaling of the 
vertical axis by the factor (1 — t). The linear relations 
between 1/ R^ and the critical temperatures are artifacts 
of GL theory, but we expect that the phase diagrams in 
FigU] are still qualitatively valid. 

In plotting FigU] P = 0.25 was used. The second or- 
der phase transition lines are independent of this ratio. 
[2^ For larger (smaller) /3, the first order transition tem- 
perature (and thus the corresponding l/iJj'**) is higher 
(lower), thus increasing (decreasing) the stability region 
for the n = phase. We also note that though our spe- 
cific calculations are for a circular disc, the n = and 
n = 1 states are distinguishable by time-reversal sym- 
metry, and so the phase transition between them should 
exist even for other geometries. 

To get more understanding of the phase diagram, we 
consider the order parameters in some detail. An exam- 
ple is as shown in Figl2] Here we have chosen 5k = 0. 
(We shall mostly be considering < (5k < 0.44 for the 
rest of this paper) . At temperatures below the first order 
transition temperature ij''* (= 0.32 here), it is a time- 
reversal-symmetry-broken state with n — Q (or n = 2). 
The flat parts of the order parameters around the cen- 
ter reflects the characteristics of the bulk system, with 
C+ 7^ but C- ~ 0. This is the preferred configuration 
at low temperatures, or equivalently, for large samples. 
For t > ij"*, the ground state becomes n = 1. The vor- 
tex structure at the center is shown in Figl^^c). The 
boundary condition ^ admits only the parallel compo- 
nent ff ■ (f) near the edge of the sample, and, being at a 
higher temperature, the radial component V-f has not 
nucleated. This implies that we have n = 1 [29|, which is 
then the preferred configuration at higher temperatures 
or intermediate size grains. The phase diagram refiects 
the competition between the boundary effect, which fa- 
vors n = I (for 5k < 0.44), versus the bulk, which prefers 
n = (or n = 2). 

At zero field, the n = (n = 2) state has surface cur- 
rent Iq in -I- (— ) (f) direction (Figl^fb)), hence a magnetic 
moment Ai along +z (— z). The n — 1 state, being time- 
reversal symmetric, has no surface current even though 
there are vortices at the center. 

Now we consider an external magnetic field B along 
the +z direction. First we consider very small fields. The 
degeneracies between n — and n — 2 are lifted, n — 
is favored since its magnetic moment A4 is parallel to the 
external field. The reduction of current near the edge due 
to external field in Figl^Jb) is due to the Meissner effect. 
For the time-reversal-symmetric state, with the applied 
magnetic field, the system has negative current near the 
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FIG. 3: b-t phase diagram: R = 1.4 and 1, Sk = 0.15,/? = 
0.25. Inset: order parameter and current at 6 = 2.5 and 2.9, 
t = 0. 
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FIG. 4: (a-c) Phase diagram: effect of rough boundary on 
the b-t phase diagram with R = 4 and Sk — 0.15. (d-f) Zero- 
field phase diagram witli rougli boundary for different Sk's. 
h = 2.0. 



edge and positive current around the vortex. 

For larger fields, the phase diagram can be modified. 
In this paper, we focus on discs with small radii and 
small magnetic field. (At larger grains or fields, cylin- 
drical symmetry can be broken due to the possibility of 
vortex lattice. We ignore this possibility in this paper). 
An example is shown for R = 1.4 in Fig|31[a). The n = 1 
state is always suppressed by the external field due to 
the kinetic energy of the Meissner current. However, for 
n = 0, due to its spontaneous magnetic moment at zero 
field, it is first enhanced by the field, then eventually 
suppressed at larger fields. Hence, above a certain field, 
n = can become more favorable than n = 1. The 
transition to the superconducting state can thus become 
reentrant as a function of magnetic field. At still higher 
fields, other n's (such as n = —1 here) can become the 
ground state. For even smaller grains, the reentry of 
superconductivity becomes more significant. An exam- 
ple is shown in Fig|31[b), where superconductivity disap- 
pears completely for intermediate fields. This reentrance 
of superconductivity is similar to the Little-Parks effect 
in s-wave superconductors in multi-connected geometries 
|30| . We note that the enhancement or suppression of 
superconductivity by B can be understood from the net 
magnetic moment of the grain, (hence the current, see 



inset of Fig 131b)) since 



dF 
dB 



-M. 



Recently, ^] studies a p-wave superconductor in a 
disc, solving the Bogoliubov-deGennes equation together 
with a weak-coupling gap equation, assuming a cylindri- 



cally symmetric Fermi surface and also a smooth bound- 
ary at R. Our results for Sk = should then be ap- 
plicable, but are different in many ways from theirs. In 
zero field and some R, they showed a transition from the 
normal state to the n — state. We however found that 
the transition from the normal state should always be 
first to the time-reversal symmetric n = 1 state, though 
the system can make a first order transition later to the 
n = state at a lower temperature for grains that are 
not too small (FigH]). While both they and we found 
reentrance of superconductivity as function of magnetic 
field, our reentrance is always to a state which is con- 
nected with the one with broken time-reversal symmetry 
(n = 0, —1, etc here), but theirs is into a state that is 
connected with the time- reversal symmetric one (n = 1) 
in zero field (Figj3]). Also, we did not find any reentrant 
behavior as a function of temperature, in contrast to j20|. 

In the above, we have assumed a smooth boundary at 
r — R. In order to mimic an imperfect boundary, we 
introduce the diffusive boundary term, Fg — KiC^/lh, 
for the free energy. |23| Here Cd = (C+ — C_)/2 and 
positive lb is the extrapolation length due to boundary 
scattering. This energy term will modify the boundary 
condition, eq.®. The factor in front of the order pa- 
rameter of right hand side becomes {K^/R) — {Ki/lf,), 
as obtained in [23 using a more microscopic considera- 
tion. Fig|3l^a-c) show the effect of rough boundary for 
R — 4. The Tc for n = 1 is suppressed more than that 
for n = 0. However, the T^'^* for the 1st order phase 
transition becomes larger for rough boundary. The sys- 
tem can be just in the superconducting state with n = 
for sufficiently rough boundary (as shown in FigHJ^c)). 
FigHfd-f) are the R-t phase diagrams for different 6k's 
without external field. Comparing these phase diagrams 
with the insets in Fig[Tl the second order transition lines 
have negative curvatures. The critical radius for the exis- 
tence of superconductivity at (5k = becomes finite. The 
critical radii at zero temperature gives a phase diagram 
similar to Fig HI except some shift of the curves to the 
left. 

In conclusion, we study the phase transition of a two- 
component superconductor in a confined geometry. We 
find that, in a large order parameter space, the system 
would exhibit multiple phase transitions. While we have 
mostly been referring to a p-wave superconducting or- 
der parameter, we expect that many of the qualitative 
features here would remain so long as the surface affects 
the two components of the order parameter differently. 
These phase transitions can be detected by, for exam- 
ple, measuring the density of states via tunneling [l6| , in 
grains of size of order of coherence length. 

This work is supported by the National Science Council 
of Taiwan under grant number NSC-98-2112-M-001 -019 
-MY3. 
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In the text, we have found that, within GL theory 
with parameters appropriate to that of a weak-coupUng 
isotropic p-wave superconductor, the transition tempera- 
ture for our circular disc with smooth boundary is not at 
all affected by the finite radius R and in fact is the same 
as that of the bulk. We shall show in this Supplemental 
Material that this result is true also within the quasi- 
classical approximation, and hence this result is not an 
artifact of the GL approximation. 

The quasiclassical green's function (4 x 4) matrix 
g{p,r,en), a function of the momentum direction p, po- 
sition r and Matsubara frequency e„, obeys 



ie„T3 - A(p,f),g +ivfp-\7g ^0 



(1) 



where we adopt the same notations as 'll|. The normal 
state green's function is given hy g — — i7rsgn(e„)r3 and 
is diagonal. A only has matrix elements that are off- 
diagonal in particle-hole space. The transition tempera- 
ture is obtained by the linearized gap equation, and hence 
it is sufficient to solve for the off-diagonal part of g to first 
order in A. Assuming spin rotationally invariant pairing 
interaction, both these off-diagonal parts of g and A have 
the same spin structure. Removing this part, we then 
simply obtain equations relating the two scalar functions 
/ and A which describe separately the off-diagonal parts 
of g and A. Eq ([T]) then reads 



(2ie„ + ivfp ■ V)/ = 2i7r(sgne„)A 



(2) 



With pairing interaction written as Vip-p' , the linearized 
gap equation reads 

A(p, r) = N[Q)TVi J2 < (p ■ p')fip', r, e„) > (3) 



where the angular bracket denotes angular average over 
p' and N{0) is the density of states at the Fermi level. 
For the bulk, A is position independent, and hence 



fiP) = TI" 



A(p) 



Eq ^ becomes 



A{p) = TTNiO)TVi XI < (p ■ 



./^A(p') 



> 



(4) 



(5) 



which determines the bulk transition temperature T°. 

For the circular disc with a smooth, spin inactive 
boundary, g obeys the boundary condition g{pin, r,en) — 
giPout, r, e„) for r at the surface, where pin and Pout de- 
note the incoming and outgoing momentum directions. 
Using the cylindrical coordinates (r, </>) for r and denoting 
the angle between p and x by 0, this boundary condition 
reads 



f{e, R, (f>, e„) = /{tt + 20-0, R, 0, e„) 



(6) 



In the text, we have seen that the highest T^ corresponds 
to the case n = 1 and the order parameters at Tc have 
the form C_^ '(r) — — C_ (r) — Cr where C is a propor- 
tionality constant. The corresponding A then has the 
form A(p,r,(/)) = fj ■ p = §r (e*^'^-^) _ e-^('^-«)) hence 

/\{p,r,(j)) = iCrsm{(l)-e) (7) 

We can verify easily that eq ^ can be solved by 



f{e,RA.^n) 



-iCrsi.n{(f) — 9) 



since {p ■ V)/ = ( cos( 
Thus we have 



a 1 sin(e- 


-0) 


d 


dr ' 




dcf, 


A{p,r) 







(8) 



/ = 0. 



/(p, f, e„) = TT 



and with this, the gap equation ^ becomes 



(9) 



A(p,f) = ',N{Q)TV^Y.< iP-p')^r^ > (^^) 



Comparison with the bulk gap equation ([5]) shows that 
this is satisfied with transition temperature Tc — T^ at 
every position r. This completes the proof. 

Thus we have verified that the transition tempera- 
ture is unchanged by the finite radius R. Since we now 
only rely on the quasiclassical approximation, this result 
should hold even for R much smaller than the zero tem- 
perature coherence length, so long as the product of the 
fermi momentum pp with R is much larger than one (and 
Tc much smaller than fermi energy so that the quasiclas- 
sical approximation holds). 
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